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We examine energy transport in an ensemble of closed quantum systems driven by stochastic 
perturbations. One can show that the probability and energy fluxes can be described in terms of 
quantum advection modes (QAM) associated with the off-diagonal elements of the density matrix. 
These QAM play the role of Landauer channels in a system with discrete energy spectrum and 
the eigenfunctions that cannot be described as plane waves. In order to determine the type of 
correlations that exist between the direction and magnitudes of each QAM and the average direction 
of energy and probability fluxes we have numerically solved the time-dependent Schrodinger equation 
describing a single particle trapped in a parabolic potential well which is perturbed by stochastic 
ripples. The ripples serve as a localized energy source and are offset to one side of the potential well. 
As the result a non-zero net energy flux flows from one part of the potential well to another across 
the symmetry center of the potential. We flnd that some modes exhibit positive correlation with the 
direction of the energy flow. Other modes, that carry a smaller energy per unit of the probability 
flux, anticorrelate with the energy flow and thus provide a backflow of the probability. The overall 
picture of energy transport that emerges from our results is very different from the conventional one 
based on a system with continuous energy spectrum. 
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I. INTRODUCTION 

Understanding the properties of the energy trans- 
port by quantum channels is one of the most important 
problems of non-equilibrium thermodynamics, important 
from the fundamental point of view as well as for its ap- 
plications. A typical problem addressed in literature is 
the energy transport between two thermal reservoirs cou- 
pled by a quantum system (a mediator). The descrip- 
tion of energy flux is usually stated in terms of Landauer 
channels described as plane waves that carry energy 
and probability flux in a given direction with occupation 
numbers determined by the temperature of the reservoir 
from which a given channel originates. 

In this paper we address the question of energy flux 
from a different perspective. One can show that both 
probability and energy fluxes can be described as super- 
positions of quantum advection modes (QAM). These 
modes are associated with the off-diagonal elements of 
the density matrix appropriately defined for our ensem- 
ble. Each such a mode carries a certain probability flux 
and a fixed, quantized amount of energy per unit of the 
probability flux. In our opinion, these modes are the pro- 
totypes of what is commonly known as Landauer trans- 
port channels. 

In order to determine the statistical properties of QAM 
and how they contribute to the overall energy and prob- 
ability flux we have considered a model of a particle 
confined within a potential well and described by the 
Schrodinger equation. We introduce a stochastic time- 
dependent perturbation of the confining potential ("rip- 
ples") that has certain well defined statistical properties. 
There is no feedback from the confined particle to the 
source of the ripples, so that the net result of the per- 
turbation is a steady increase in energy of the confined 



particle. 

Since the power source is localized, there is a nonzero 
average energy flow from one spatial region of the con- 
fining potential to another. Thus, we can study the sta- 
tistical properties of the energy and probability fluxes 
without addressing a more involved problem of two ther- 
mal reservoirs entangled with the mediator. This model 
is somewhat similar to that describing trapped ions and 
atoms where electromagnetic perturbations lead to grad- 
ual heating of the confined particles [1-0| • 

From the solution vE'(r, t) of the Schrodinger equation, 
which includes perturbation by ripples, one can calculate 
the probability flux jp (r, t) . Similarly, one can introduce 
the energy density t(f,t) and the energy flux jE{r,t), 
both expressed in terms of vE'(r, t) and its time derivative 
^'(r, t), that satisfy the energy conservation condition in 
the form of the continuity equation. This equation in- 
cludes also the energy source determined by the time- 
dependent perturbation. The source of stochasticity in 
our approach is a classical Brownian particle immersed in 
a thermal bath which is characterized by temperature T. 
The motion of the Brownian particle creates ripples in the 
confining potential. The spatial trajectory of the Brow- 
nian particle is determined by the standard Langevin 
equation. Thus, the quantum particle described by the 
Schrodinger equation is affected by the movement of the 
classical Brownian particle without any feedback. 

Our approach is as follows: We generate, solving nu- 
merically the Langevin equation, a time dependent tra- 
jectory of the classical Brownian particle. Then, feeding 
this trajectory into the Schrodinger equation as part of 
the perturbation potential, we obtain the solution 4'(r, t). 
After that we start the procedure all over again, gener- 
ating another random trajectory and another solution 
^'(r, t). By repeating this process many times we ac- 
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cumulate an ensemble consisting of the wave functions 
'i>{f,t). Each of these solutions has evolved under the 
action of the perturbations, all of which have the same 
statistical properties determined by the movement of the 
Brownian particle. It is important to emphasize that this 
ensemble is not a traditional thermodynamic ensemble, 
it is rather a collection of evolving pure quantum states. 

For each member of this ensemble we calculate proba- 
bility and energy fluxes and analyze the statistical prop- 
erties of the resulting ensemble of fluxes. 

This paper is structured as follows. Section II defines 
the quantum advection modes in terms of the off-diagonal 
elements of the density matrix. Section III introduces the 
model which we use to determine the statistical proper- 
ties of these modes. Section IV presents the results which 
include the distributions of the probability and energy 
fluxes. The ensemble average energy flux is treated as a 
superposition of two statistically orthogonal fluxes. One 
is a purely advective energy flux fully correlated with the 
probability flux. The other component - the thermal flux 
- has nonzero average value and corresponds to the en- 
ergy flow uncorrelated with the movement of the particle. 
The subsection IV(B) presents the correlations between 
the direction of the individual QAM and the direction of 
the energy flow averaged over the time intervals where 
the net probability flux is zero. 



II. QUANTUM ADVECTION MODES 

Let us first consider some general statements regarding 
energy fiux that corresponds to a wave function 4'(r, t) 
satisfying a time-dependent Schrodinger equation: 

where p = ~ihW and the potential V{f,t) = Vo(r) + 
U {r, t) consists of a confining potential Vq (r ) and a 
superimposed upon it time dependent perturbation U. 
The probability of finding a particle in a given place 
P{r,t) — |^'(r, i)p obeys the continuity equation 



dP 
'di 



+ V • = 0, 



(2) 



where the probability fiux jp is given by 

= los{^*v^). (3) 

TO 

The energy density can be defined as follows 
e(r,i) = ]^[^*{H^) + (i?**)*] = ^{^*H^). (4) 

It is a real number, which after integration over the total 
volume gives the expectation value of the Hamiltonian 



Etot = / €{r,t)dr^ {H). 



(5) 



By virtue of Eq. (1) the energy density can be written 
in the form 



(6) 



The quantity defined by Eq. (6) was introduced by D. 
Bohm [Sj . He suggested to readers to prove that it satis- 
fies the conservation law (in the case of time-independent 
Hamiltonian). At the risk of upsetting future students 
by denying them a pleasure to derive the energy continu- 
ity equation themselves, we give in the Appendix a brief 
derivation of the following equation 



de 
di 



^■h^Q{r,t). 



(7) 



Here je is the energy fiux 

h 



2to L J 
A more convenient for numerical calculations form is 



jE 



2m I 



(8) 



(9) 



The energy source/sink in Eq. (7) is determined by the 
time dependent part of the Hamiltonian, 

Q = T>|*|2. (10) 

Any solution of Eq. (1) can be presented in a form 



-tpnir); iOn=En/h. (11) 



Here we consider the wavcfunctions ij^n (?^) to be the eigen- 
states of the unperturbed Hamiltonian 

Ho = ^ + Vo{r)- Ho^n = E„4,r^- (12) 

Substituting (11) in Eq. (3) we obtain the probability 
fiux expressed in terms of the stochastic coefficients c„ (t) 



jp{f,t) = — V 

TO ^-^ 



*) 5„„(f), (13) 



where LUnm = — '^m > 0. The spatial distribution of 
the flux is determined by the generalized Wronskians of 
the two wave functions 



(14) 



This expression for the probability fiux is universal 
even though the wave functions tpni^ are not the eigen- 
states of the full Hamiltonian H. It is a bit more com- 
plicated for the energy fiux. Let us consider a situation 
when the time dependent perturbation U (r, t) is local- 
ized, so that there are spatial regions within the confining 
potential free of perturbation. In these regions 



(15) 
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and, correspondingly, 

H{f,t)Mr^ ^ EnMr^- (16) 

Within these free of perturbation regions we obtain from 
Eqs. (8) and (11) 

(17) 

Thus, 

Jp{r,t) ^ qnni{f,t); (18) 

jE{r,t) = ^ enmqnm{r,t). 

Each of the quantum advection modes (QAM) 

qnm{r, t) = {cnc*^e-"^-"^') gnm{r^ (19) 

carries a certain amount of the probability flux deter- 
mined by the correlator (c„c*je~'""'"*). Wc will loosely 
call these correlators off-diagonal elements of the density 
matrix, even though the ensemble over which the aver- 
aging takes place is not necessarily a conventional ther- 
modynamic ensemble. The amount of energy carried by 
each mode per unit of the probability flux is fixed, given 
by the average energy of the two levels involved. 



Hereafter we will ignore the spin degeneracy. We 
will call two modes qnm and qn'm' degenerate if they 
carry the same energy per unit of the probability flux, 
namely if e„m = Cn'm'- Consider, as an example, 
one-dimensional oscillator with equidistant energy lev- 
els En = hjo{n + 1/2), so that e„m = fiw(r?, + m + l)/2. 
The lowest energy mode with positive spatial parity qn) 
is nondegcncratc and carries one quantiim of energy per 
unit of the probability flux. The lowest energy mode with 
negative parity q2o carries fractional amount itvjj/2 per 
unit of the probability flux and is also nondegenerate. 
Two degenerate modes, both of positive parity, q^Q and 
(721 carry two quanta of energy per unit of the probability 
flux, etc. 

A few immediately obvious conclusions follow from 
this analysis. Each QAM corresponds to purely advec- 
tive energy transfer, namely the energy flux carried by 
each mode is directly proportional to the probability flux. 
However, in a system where the net average probability 
flux is zero (there is no net mass or charge transfer) it is 
possible to have a non-zero energy flow, provided that at 
least two non-degenerate QAM are activated and both 
have nonzero average, so that jp = qa + qp = 0, but 
js = ^aqa + ^0q0 7^ 0. Two nondegenerate modes re- 
quire at least three energy levels involved. A two-level 
system cannot conduct energy without carrying nonzero 



probability (mass or charge ) flux. The energy flow un- 
accompanied by mass flow is usually defined as thermal 
conduction (heat flux). 

It is important to realize also that energy flow through 
a quantum system (such as a qubit) establishes a degree 
of persistent coherence between several cncirgy levels, so 
that at least some correlators in (19) do not average to 
zero for as long as the external conditions that create the 
energy flow are maintained. 



III. MODEL 

In order to determine the properties of the quantum 
advection modes we have studied a following model. We 
have solve numerically the one-dimensional Schrodinger 
equation (1) with parabolic conflnement potential 

Vo{x) = . (21) 

A time dependent perturbation is introduced as a fixed 
shape "protrusion", U{x — xi{t)), riding on the back of 
a classical Brownian particle whose trajectory is Xi{t), 
Fig. 1(a). Specifically, we choose 

Uix - x,it)) = C/oexp . (22) 

The classical Brownian particle is trapped in its own har- 
monic potential well and is in equilibrium with the ther- 
mal bath of temperature T. It is important to emphasize 
that the center of the potential well in which the Brown- 
ian particle is trapped is located off center of the potential 
well Vo(a:), as shown in Fig. 1(a). We can choose parame- 
ters, such as the location of the center of the trap and the 
stiffness constant of the trap in such a way that most of 
the time the Brownian particle and, correspondingly, the 
perturbation are localized on the right-hand side, x > 0, 
of the confining potential Vq{x). Thus, the left-hand side 
of Vo{x) is mostly free of perturbation, which is a con- 
dition, Eqs. (15) and (16), necessary for defining the 
energy flux in terms of the quantum advection modes, 
Eqs. (17) and (18). 

To make the perturbation a ripple, rather than a 
protrusion, we define it by subtracting the static part 
U{x — a). Thus, the potential in Eq. (1) is defined as 

V{x,t) = Va{x)+SU{x,t); (23) 
6Uix,t) = U{x - xi{t)) - U{x - a). 

Here a > is the location of the center of the trap for 
the Brownian particle. If the Brownian particle were at 

rest at the bottom of its trap, SU = 0. Since its motion 
is determined by the temperature of the thermal reser- 
voir, the perturbation vanishes at T = 0. At all times 
jTao (^' *)'^'^ ^ ^- Figui"c 1(b) shows a snapshot of the 
total potential Vo{x) + 5U{x,t). 
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FIG. 1. (a) A sketch of the potential V{x,t) consisting 
of a static, symmetric confining potential with super- 
imposed protrusion driven by a classical Brownian 
particle immersed in a thermal bath, (b) A snapshot 
of the potential V{x,t), Eq. (23). 



Thus, we solve numerically the folowing Schrodinger 
equation 



dt 2m dx^ 



+ {Vo{x)+SU{x,t))'i>. (24) 



A trajectory of the Brownian particle necessary to 
define the perturbation (23) is obtained by solving nu- 
merically the Langevin equation which determines the 
stochastic trajectories xi{t), where 

Mxi+Tjxi+Mujfixi-a) ^ F, (25) 

with the random force correlator 



{F{t)F{t')) ^2rjkBTS{t-t'). 



(26) 



Our approach to analyzing the processes caused by 
such stochastic perturbation is as follows. A trajectory 
xi{t) is generated by solving numerically Eqs. (25, 26). 
This trajectory is then fed into the numerical algorithm 
which solves Eq. (24) with the potential given by (22, 
23). This gives us the wave function ^'i(x,t). Then an- 
other trajectory xi{t) is generated by the same equations 
(25, 26) and the second solution of Eq. (24) ^'2(0;,^) is ob- 
tained. This protocol is repeated TV times. As the result 
we obtain an ensemble of solutions {^i(x, t) . . . '^'^{x, t)} 



which allows us to determine the ensemble- and time- 
average, as well as the distributions, of different physical 
quantities such as the probability and energy fluxes, the 
density matrix elements, etc. The data presented here are 
the result of averaging over iV = 10^ copies of the quan- 
tum system exposed to the action of the stochastic rip- 
ples. It should be noted that this is an ensemble of pure 
quantum states, so that it contains all information about 
the quantum system, rather than a truncated amount of 
information described by the conventional density ma- 
trix. The total energy of this ensemble is not fixed, in 
fact it increases at a certain rate. In order to analyze our 
results we do not introduce any "average" wave func- 
tion. We are only averaging the probabilities and fluxes, 
i.e. the quadratic forms of the individual wave functions 
that constitute this ensemble. 



A. Implementation 

In the appropriate dimensionless variables 
T — Lot; ^ — x/£, 

Eq. (24) takes form 



where 



£ = 



1/2 



(27) 



(28) 



(29) 



Energies and temperature are measured in their natural 
unit hui, so that 



huj 



The unperturbed energy levels are given by 

En=n+ 1/2, 
and dimensionless temperature is defined as 

The Langevin equation (25) takes form 



(30) 



(31) 



(32) 



(33) 



where ^1 = (xi — a) / 1. The relative frequency of the con- 
fining potential and the relative damping rate are given 
by 



7] 



2 _ 



(34) 
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The third parameter describing the Brownian particle is 
the relative stiffness of its trap 



k muj 



ki Mujf 



(35) 



which determines the characteristic scale of the displace- 
ment of the Brownian particle from the center of its trap 



2 fcsT_ il 



-; -p; — Kb 



fci ' £2 

The random force in Eq. (33) is defined as 

F 



with the correlator 



(/(T)/(r')) - 2j0{Kn')S{T 



(36) 



(37) 



(38) 



The Langevin equation (25, 26) describes a particle in 
thermal equilibrimxi. Therefore, the trajectories obtained 
as its solutions must have the statistical properties de- 
fined by the Boltzmann probability distributions for the 
displacement and velocity 



P(xi) 



2TTkBT 



1/2 



exp 



MijjI{xi - af 
2fc^T 



and 



P(ii) 



M 



27rkBT 



1/2 



exp 



2kBT 



(39) 



(40) 



Here P{xi)dxi and P{xi)dxi are relative frequencies 
with which the trajectories acquire the respective values. 
We have used these conditions as a way to validate the 
correctness of our numerical solutions of the Langevin 
equation. 

The numerical parameters that were used in this study 
were as follows. The temperature range was ~ 1. The 
frequency of the trap uji was chosen close to the main 
frequency cj, and the damping rate 7 was of the same 
order of magnitude as ft. The data set shown below 
correspond to 



n = 1.1; 7 = 1.6. 



(41) 



If we choose the range of parameters substantially dif- 
ferent from that shown above, the perturbation will be 
either too weak or too strong to obtain reliable results 
using numerical methods of solution of the Schrodinger 
equation. 

We also need to ensure that the probability of the per- 
turbation (23) to extend to the left-hand side {x < 0) 
of the confining potential Vo{x) is statistically negligible. 
This is necessary to meet conditions (15, 16) so that we 
can use the formalism of QAM to describe the transport 
in the left-hand side of the potential well. For that pur- 
pose the parameters of the trap displacement a and the 



stiffness of the trap have to be chosen so that the prob- 
ability to find the Brownian particle at xi < was neg- 
ligibly small. The probability (39) to find the Brownian 
particle at xi = is given by 



/ 1 



1x1=0 



V 27rK6' 



1/2 



exp{-a^/2fKe}. (42) 



In our numerical calculations we have used the values 



a/£= I; K= 1/36. 



(43) 



If in Eq. (23) the displacement a/£ is substantially 
greater than 1, it reduces exponentially the coupling con- 
stants between the energy levels. The stiffness of the trap 
36 times greater than that of the potential Vq, coupled 
with the width of the Gaussian (22) 



6 = 1/8 



(44) 



ensures that the left-hand side of the potential well Vq (x) 
is practically free of perturbation. 



B. Numerical methods 

To solve the Langevin equation given by Eq. (33), we 
used the standard fourth-order Runge-Kutta (RK4) iter- 
ative method. The random force, Eq. (38), was defined 
as follows: 



Ati/2 ' 



(45) 



where At is the time step and 77 is a pseudo-random vari- 
able distributed with the probability 



1 



1/2 



exp{~r]^/4-f9Kn^}. (46) 



Validation of the algorithm was carried out by comparing 
the statistical properties of the trajectories with condi- 
tions (39, 40). 

The solution for time-dependent Schrodinger equation 
given by Eq.(28) is obtained using the Crank- Nicolson 
method, similar to that discussed in [9l-[ll|. In each in- 
dividual run the perturbation was turned on and turned 
off gradually at the beginning and the end of the run in 
order to avoid the excitation of too many energy levels. 
Validation of the algorithm was carried out by control- 
ling the norm of the wave function and by running the 
algorithm without perturbation and comparing the re- 
sults with the exact solution. It was determined that 
error greater than 5% accumulates after approximately 
20 — 30 classical periods 27r/w. This places the limit on 
the duration of time that we have used to run numerical 
calculation for the individual solution 5'(x,i). 
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FIG. 2. (a) Ensemble average energy, Eq. (47), for 
three different values of the friction coefficient 7 in 
Eq. (33). (b) Two examples of the energy of an indi- 
vidual member of the ensemble Elot • In contrast to 
the ensemble average energy, the energy of an indi- 
vidual quantum system can decrease over some time 
intervals. 



IV. RESULTS 



Figure 2 (a) shows the ensemble average energy 



N 



(47) 



1=1 



obtained by averaging over N = lO"* sohitions ^i{x,t) of 
Eq. (24). The individual energies El^^ were calculated 
according to Eqs. (5, 6) 

/oc 
^{^*{x,t)i',{x,t))dx. (48) 
-00 

All results shown in this paper correspond to the 
ground state initial condition for all members of the en- 
semble, so that Elgf{0)/haj — 1/2. On average, the per- 
turbation continuously pumps energy into the system. 
This is true in the case of a spatially uniform perturba- 
tions (heating of trapped ions) 0, hj and in our model 
as well. The rate of energy increase depends on the pa- 
rameters of the Brownian particle. Other things being 




0.25 



FIG. 3. (a) Time dependence of the ensemble average 
Ei/2 and P1/2. 



equal, the lower friction leads to greater rate of energy 
increase. However, the averaging conceals the full pic- 
ture of this phenomenon. Figure 2 (b) shows the energy 
of two arbitrary chosen individual members of the ensem- 
ble. On the level of an individual quantum system the 
energy does not monotonically increase. There are inter- 
vals of time when the perturbation serves as the energy 
drain. 



A. Probability and Energy Fluxes 

In our model there are no two heat reservoirs between 
which energy can flow. Nevertheless, it is easy to see 
that this model does allow a persistent energy flow in 
the absense of the net probability flux. Let us consider a 
"half-bucket" energy 



El 



/2 



e{x, t)dx. 



If we integrate Eq. (7), we get 
dEi/2 



dt 



= -JE{t) 



(49) 



(50) 



Here we have taken into account that the energy source 
Q{x, t) is localized at x > 0, so that the energy inflow into 
the X < half of the confinement is due to the energy 
flux across the centerline. In other words, the energy flux 
can be calculated on the basis of Eq. (9) or, in the spatial 
regions free of perturbation, on the basis of Eq. (50). 

Figure 3 shows the ensemble average i?i/2. The av- 
eraging procedure is the same as in (47). Obviously, 
the excess of energy introduced by the perturbation does 
not accumulate in one half of the confining potential, so 
that -Ei/2 ~ £'tot/2 almost monotonically increases and, 
therefore, there is a persistent net energy flux across the 
centerline (from right to left) Je < . 

Similarly we can consider 



Pl/2{t) 



P{x, t)dx, 



(51) 
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FIG. 4. (a) The histogram shows the distribution of the values of the energy flux between two halves of the 
confining potential, at a:: = 0. (b) Time dependence of the ensemble average of the energy flux. Notice the scale 
difference between the values of je in Fig. 4(a) and 4(b). (c) The histogram of the distribution of the values 
of the advective component -qip of the energy flux between two halves of the confining potential, (d) Time 
dependence of the ensemble average of the advective component of the energy flux, (e) The histogram of the 
distribution of the values of the thermal component jth of the energy flux between two halves of the confining 
potential, (f) Time dependence of the ensemble average of the thermal component of the energy flux. 



which determines the probabihty to find a particle at 
X < 0. Integrating Eq. (2) we get 



dPi 



/2 



dt 



-hit) 



(52) 



x=0 



Figure 3 also shows the time dependence of ensemble av- 
erage Pi/2- Obviously, the probability does not accumu- 
late in one half of the confining potential, so that there 
are numerous time intervals over which the time aver- 
age probability flux is zero, (^1/2(^2) — A/2(ii))/(i2 — 
ti) = 0, but the time average energy flux is not zero, 
(i?i/2(t2) — £'1/2(^1 ))/(i2 - ti) 7^ 0. One can notice that 



there is small imbalance between the two halves of the 
potential well. The region where the perturbation is lo- 
cated becomes somewhat "deeper", so that on average 
Pi/2 < 0.5, but the imbalance is small 



Pi 



1/2 



0.5 



0.5 



-0.02. 



(53) 



For the ensemble of solutions 
{'i'i{x,t) . . .'i>N{x,t)} we calculate 

ing ensembles of fluxes {jp^\x,t) 
,■(1) 



.j'j^\x,t)} using Eqs. 



of Eq. (24) 
the correspond- 
.jl^\x,t)} and 
(3) and (9) at any 
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given point as functions of time. 

A histogram in Fig. 4 (a) shows the normaUzcd to 
unity probabihty distribution dW{jE)/djE for the energy 
flux across the center line x = 0. The natural units for 
the energy flux in one-dimensional system is hu)'^. This 
distribution is not a normal Gaussian bcicause the tails 
are "fatter" than normal. It can be described as a super- 
position of normal and Laplace distribution 



dW{x) - cxp(-|x|). 



(54) 



Figure 4(b) shows the ensemble average of the energy 
flux across the centerline. The averaging procedure is 
the same as in (47) 



N 



i^(i) = Ar-i^j«(t) 



(55) 



x=0 



As expected, one can sec that there is a negative bias in 
this stocliastic quantity corresponding to the net energy 
flow from the right-hand side of the potential well where 
the perturbation is localized to the free of perturbation 
left-hand side. Notice the difference in scale between Fig. 
4(a) and 4(b). The ensemble average of the energy flux 
is relatively small 



4^ -±(0.02-0.04), 



(56) 



while the energy flux in an individual member of the en- 
emble, as evident from Fig. 4(a), is substantially greater 



Je 



±0.5. 



(57) 



The large stochastic fluctuations of Je within the ensem- 
ble cancel each other out, resulting in relatively small 
ensemble average net energy flow. This is the reason the 
distribution in Fig. 4(a) looks almost symmetric. 

It is important to identify two different types of energy 
transfer. The center of mass of the particle, as deter- 
mined by the probability flux, randomly moves between 
two halves of the confining potential and carries with it 
a certain amount of energy. The average of this com- 
ponent of the energy flux is obviously zero. The energy 
flux of a different nature, uncorrelated with the particle 
movement, we will call, for lack of a better term, the 
thermal component of the energy flux. For each member 
of the ensemble the energy flux can be presented as a 
superposition of two stochastic components 



Je = Tjjp +jth, 



(58) 



where rjjp is the advective component of the energy flux. 
In a steady state ensemble, r] would be a constant - the 
average energy of advective flow. In our ensemble with 
average energy rising, rj slowly increases with time. The 
second component jth is the "thermal flux". We define 
the thermal flux as a part of the energy flux Je that is 
statistically orthogonal to the probability flux, namely 



{{jpjth)) = 0. 



(59) 



Here we use double brackets to indicate a stistical aver- 
aging, which may include both time and ensemble aver- 
aging, and differentiate it from the quantum average as 
in Eq. (5). The condition of orthogonality is 



V = 



{{jpjE)) 



Because of the orthogonality only one of the components, 
either the advection or the thermal component of the to- 
tal energy ffux, can have a non-zero average. In onv case 
the average probability flux is zero, so that the net en- 
ergy flux is part of the thermal component jth- The ad- 
vective component describes the energy flow associated 
with random sloshing of the particle around the poten- 
tial well and is not associated with the heating of the 
system. It should also be noted that the condition of 
orthogonality, Eq. (59), is equivalent to the condition of 
minimization of the difference je — rjjpi namely 



mpf)) = 0- 



(61) 



In Fig. 4(c) the distribution of the; advective compo- 
nent is shown. This is a purely symmetric distribution. 
The value of the average energy 77 carried by the probabil- 
ity flux increses with time as the total energy increases. 
Figure 4(d) shows the time dependence of the ensemble 
average advective component of the energy flux. 

Figure 4(e) shows the distribution of the thermal com- 
ponent. It reminds more closely an asymmetric Laplace 
distribution 



dW{j 



th 



djth 



■ exp(-|j(/j|/j±) 



(62) 



with the the constant < j_ . Even in the thermal com- 
ponent of the energy flux the fluctuations are substan- 
tially greater than the net average, so that \j+ — j-\ <SC 
{j+ + j_)/2. In Fig. 4(f) the ensemble average of the 
thermal component (obtained similar to (55)) is shown. 
Here the fluctuations cancel each other out and the neg- 
ative bias is very clearly pronounced. 

In our model the amplitudes of fluctuations of the en- 
ergy flux arc substantially greater than its net average. 
Segregation of the thermal flux component from the to- 
tal energy flux allows us to better resolve a small effect 
(non-zero net energy flux) on the background of the large 
fluctuations. 



B. Quantum Advection Modes as Energy Carriers 

One of the main goals of this study is to determine how 
different QAM contribute to the overall energy flux. As 
one can expect, some of the modes will carry probability 
(charge and/or mass) and energy in the direction of the 
average energy flux, while the other modes will provide 
the backflow and carry probability in the opposite direc- 
tion. Thus, all these advection modes can be identified 
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FIG. 5. (a) Time dependence of the half-bucket prob- 
ability for an individual member of the ensemble, Eq. 
(51). (b) Time dependence of the half-bucket energy, 
Eq. (49). Four consecutive averaging time intervals 
determined by the condition (63) are shown by the 
dashed lines. These intervals are different for every 
member of the ensemble. 



as the Landauer channels in the quantum system with 
discrete energy spectrum. 

As an example, let us consider conditions when the av- 
erage, over a certain time interval, probability flux across 
a given point is zero, but the average energy flux over the 
same time interval is non-zero. Furthermore, let us con- 
sider, for the sake of argument, that there are only two 
active QAM, Qa and and < e^. Then, obviously, the 
fact that the average jp = Qa + = requires that the 
direction of the low energy mode, qa, anticorrelates with 
the direction of the energy flux js = £aqa + £i3qp 7^ 0. 
The lower energy mode provides the backflow of the 
probability. The direction of the mode qp which carries 
greater energy per unit of the probability flux must cor- 
relate with the direction of the energy flux. When there 
are several activated modes, the correlation between the 
average direction of the energy flux and each individ- 
ual mode becomes less than 100%, but it still must be 
present. 

Figure 5(a,b) shows a typical time variation of the half- 
bucket probability and energy, Eqs. (49) and (51), of an 
individual member of the ensemble. There are several 
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FIG. 6. (a) The distribution of the values of the en- 
ergy flux in the subensemble defined by a condi- 
tion < 0. (b) The distribution of the values of the 
lowest QAM qlo in the same subensemble. (c) The 
distribution of the values of the QAM ^21 in the same 
subensemble. (d) The distribution of the values of 
the QAM q-jQ in the same subensemble. 
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Table 1 





qio 


q2i 

Enm=21ia) 


qso 


q32 

£nm=31ia) 


q4i 

£nm=3TKjO 


qso 


q43 


q52 

£nm=41ia} 


q54 

£nm=51ia) 




-0.8213 


0.3875 


0.3410 


0.5134 


0.3770 


0.1542 


0.2213 


0.3618 


0.1968 


Rnm 


-0.8547 


0.4324 


0.4049 


0.5459 


0.2797 


0.1705 


0.1560 


0.2559 


0.2252 


^nm 


-0.8711 


0.4567 


0.3900 


0.5771 


0.3059 


0.1339 


0.1117 


0.2220 


0.2087 




-0.8794 


0.4365 


0.3868 


0.6246 


0.2488 


0.1508 


0.1606 


0.2317 


0.1983 



FIG. 7. Table 1: The values of the correlators Rnm, Eq. (68), for the first nine modes (jnm- The top row 
corresponds to the earliest time intervals. The other rows show the values of the correlators at progressively 
more distant times. 



instances when 



Pl/l{tk) = Pl/2{tk+l), 



(63) 



so that the average, over the time interval tk+i —tk, prob- 
ability flux is zero. However, the average energy flux over 
the same interval is, in most cases, nonzero 



1 



tk+1 — tk 



j-tk+i 
Jtk 



dtfEit) 



= (64) 



x=0 



1/2 V 



tk+1 — tk 



Here the superscript i = I . . . N indicates a member of 
the ensemble characterized by a wave function ^'i(x,<:). 
The time intervals tk,tk+i are different for different mem- 
bers of the ensemble. We have used several consecutive 
intervals progressively farther from the starting moment. 
The duration tk+i — tk of the averaging intervals is of 
the order of one classical period, oj{tk+i — tk)/2TT ~ 1. 
For each member of the ensemble we determine the time 
intervals satisfying condition (63) and calculate the time 
average energy flux given by Eq. (64) as shown in Fig. 
5. 

The average value of different QAM is calculated for 
each member of the ensemble by calculating the ampli- 
tudes Cn{t), Eq. (11), and averaging the contribution of 
each mode 



^it) 



h 



-ioJnmt 



x=0 



over the the same time interval as in (64), 

(lnm{t)dt. 



tk+1 — tk Jt 



(65) 



(66) 



At a; = only odd-even modes provide nonzero contri- 
bution because 



= V'n(a;)V'm(a;) - i;n{x)tp'm{x) 



(67) 



vanishes at x = if the indices n, m are both odd or even. 
Only the modes with positive spatial parity contribute to 
the flow of matter and energy across the symmetry point. 

A quantitative measure of the overall correlations be- 
tween the energy flux (64) and a given mode (66) is de- 
termined by the correlators 



P'nm — 



nm 



where 



JV 



N 



(Qnm) = 



01) 



(68) 



(69) 



Figure 6 provides a qualitative view of the correlations 
between the energy flux and different modes. In Fig. 6(a) 
we show a histogram of a subensemble of the values of j|; 
(Eq. (64)). This subensemble includes over half of the 
total number of members of the full ensemble for which 
Je — 0- The rest of the ensemble has > 0. Figure 
6(b) shows the distribution of the values of the lowest 
QAM qIq for the same subensemble. This mode strongly 
anticorrelates with the direction of the energy flow and 
on average is responsible for the backflow of substance 
(probability flux) and energy in the direction opposite to 
the direction of the average energy flow. We see exactly 
the same anticorrelation between and qIq in the rest 
of the ensemble for which > 0. 

Figure 6(c) shows the distribution of the values of ^21) 
one of the two degenerate modes which carry two quanta 
of energy per unit of the probability flux for the same 
subensemble. This mode clearly correlates with the en- 
ergy flux and on average carries probability and energy in 
the same direction as the total energy flow. Figure 6(d) 
shows the distribution of the second degenerate mode qIq. 
For this mode the distribution appears more symmetric, 
so that its contribution to the net energy flow is some- 
what smaller than that of qli- Notice, that for qIq and 
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^21 the frequency w„m = u), see Eq. (19), while for q^gg 
W30 = 3aj. 

Table 1 shows the values of the correlators Rnm (Eq. 
(68)) for nine lowest modes which account for about 90% 
of the total contribution to the values of j^. The av- 
eraging takes place over the whole ensemble. For each 
member of the ensemble 

j; = E ^nm = 0; (70) 

n>m 

n>m 

The rows in the table correspond to successive time in- 
tervals like the ones shown in Fig. 5. The top row corre- 
sponds to the earliest time intervals and the other rows 
give the values of the correlators at progressively later 
time. 

The values of the correlators confirm the qualitative 
picture seen in Fig. 6. The strongest correlation is ex- 
hibited by the lowest energy mode g^Q, which appears to 
be single-handedly carries the backflow of the probability 
flux equal to that carried by all other modes in the oppo- 
site direction. This is the only mode that anticorrelates 
with the direction of the net energy flow in the large ma- 
jority of the ensemble members. All other higher energy 
modes exhibit positive correlations with the direction of 
the energy flow. We expected that the lowest mode will 
be responsible for the backflow, but it is rather surprising 
that all other modes exhibit positive correlations within 
the statistical margins of error. The time progression 
does not change these correlators in any significant way. 

It is possible, that the choice of the ground state as the 
initial condition, which makes the correlator coc\ greater 
than all other correlators, makes it inevitable that gio 
is the only mode responsible for the backflow. All other 
modes need to combine their contributions constructively 
in order to cancel that of qio to the total probability flux 
and prevail in determining the direction of the energy 
flux because they carry greater energy per unit of the 
probability flux. 

V. SUMMARY AND SPECULATIONS 

The main purpose of this paper is to clarify the no- 
tion of Landauer channels in a system with discrete en- 
ergy spectrum, whose eigenfunctions cannot be reduced 
to plane waves. We also wanted to address a situation in 
which there is no net probability (charge and/or mass) 
flux, but there is a net energy flux. In our model we pro- 
duce the persistent energy flux by pumping energy into 
the system with the help of a localized and asymmetri- 
cally spaced power source. 

In the region free of the time-dependent perturbation 
the energy and probability flux can be defined as a su- 
perposition of quantum advection modes, each associated 
with an off-diagonal element of the density matrix. Our 
results show that these modes play the role of Landauer 



channels, but in a way very different from their conven- 
tional interpretation. First of all, each mode does not 
simply carry probability and energy in a certain direc- 
tion. The direction and amplitude of the modes fluctu- 
ate, and only statistically correlate with the net energy 
flow. The lowest energy mode provides on average the 
backflow of the probability and energy. The higher en- 
ergy modes have positive correlations with the energy 
flux. Thus, the total probability flux can be zero because 
the contribution of the backflow mode cancels out that 
of the higher energy modes. However, since the higher 
modes carry greater amount of energy per unit of the 
probability flux, they determine the value and the direc- 
tion of the net energy flow. 

Usually, the energy flow unaccompanied by the mass 
flow is deflned as heat transfer. The results shown above 
indicate that a constant heat flux exists due to persistent 
coherence between the amplitudes of the different energy 
levels. As long as the external conditions that induce the 
persistent energy flow are maintained, the ensemble av- 
erage values of the correlators c„c*je~*'^"'"* are not zero. 

This observation might have implications beyond the 
heat transfer topic. Decoherence of qubits, the funda- 
mental elements of the hypothetical quantum computers, 
is a major obstacle on a way to their construction. One 
might suggest that a qubit placed in contact with the 
nonequilibrium environment, such that it induces a per- 
sistent energy flow through the qubit, will be protected 
to some extend from complete decoherence. One can 
see the illustration of this point in Table 1. The mode 
qio ~ S (ciCge"*'^""*) has over 80% correlation with the 
net heat flux. Thus, strong persistent coherence between 
the ground and the first excited state exists for as long as 
we maintain the energy fiow through this quantum sys- 
tem. We should mention again that trapped ions, one of 
the early candidates for the role of qubits, are well de- 
scribed by a model of a particle confined to a harmonic 
potential well and exposed to stochastic perturbation due 
to electromagnetic noise [1-0|- Of course, at this point 
this is only a speculation and more work needs to be done 
to determine whether this hypothesis has any merit. 
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VII. APPENDIX: DERIVATION OF ENERGY 
CONTINUITY EQUATION 

There are many ways one can crumble this cookie. 
Here we will show that Eqs. (7) and (9) are correct. The 
energy density deflned by Eq. (4) and the Hamiltonian 
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(1) is given by 



2m ^ /II 



The rate of change 



where 



and 



de . 



ei = 3? ( **V2* + 4'*V^* ) 

2m V / 



e2 = y|*|2 + F 



We can rewrite ei in the form 



Then, it is easy to see that 

sfi ( ^*v2^f - **V^* I = 

2m \ / 

-— 3i|V-(**V*-*V**)| =-V-jE. 
Now we only need to show that 

dt m ^ ^ 

This condition is equivalent to the probability conserva- 
tion condition (2). Indeed, by virtue of Eq. (1) 

SR{**V^*} = — 9(iI**V^*) = 

m m 

— y3(**V2*) = — F3{V • (^'*V^')} = VV ■ t. 

mm 

If we take the wave function \1/ in the form 
*(r,i) = Re^^, 



so that 



we get 



jp = -R^VS, 
m 



fe2 . . 

ei = — -5R 1 **V^* - **V^* + 2**V^* ) 
2m V / 



3E = 7^{R^R - R^R} - hSjp. 

2m 
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